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Abstract 


This paper presents a review of high-order and optimized finite-difference methods for nu- 
merically simulating the propagation and scattering of linear waves, such as electromagnetic, 
acoustic, or elastic waves. The spatial operators reviewed include compact schemes, non- 
compact schemes, schemes on staggered grids, and schemes which are optimized to produce 
specific characteristics. The time-marching methods discussed include Runge-Kutta meth- 
ods, Adams- Bashforth methods, and the leapfrog method. In addition, the following fourth- 
order fully-discrete finite-difference methods are considered: a one-step implicit scheme with 
a three-point spatial stencil, a one-step explicit scheme with a five-point spatial stencil, and 
a two-step explicit scheme with a five-point spatial stencil. For each method studied, the 
number of grid points per wavelength required for accurate simulation of wave propagation 
over large distances is presented. Recommendations are made with respect to the suitability 
of the methods for specific problems and practical aspects of their use, such as appropriate 
Courant numbers and grid densities. Avenues for future research are suggested. 


Introduction 

Numerical simulation can play an important role in the context of engineering design and 
in improving our understanding of complex systems. The subject of wave phenomena, in- 
cluding electromagnetic, elastic, and acoustic waves, is a promising area of application for 
such simulations. Lighthill [1] and Taflove [2] discuss the prospects for computational aeroa- 
coustics and electromagnetics, respectively. The computational requirements of accurate 
simulations of the propagation and scattering of waves can be high, particularly if the size of 
the geometry under study is much larger than the wavelength. Consequently, there has been 
considerable recent effort directed towards improving the efficiency of numerical methods for 
simulating wave phenomena. 

In electromagnetics, the most popular approach to the numerical solution of the time- 
domain Maxwell equations for numerous applications has been the algorithm of Yee [3], which 
was named the finite-difference time-domain (FDTD) method by Taflove [4]. This algorithm 
combines second-order centered differences on a staggered grid in space with second-order 
staggered leapfrog time marching. Its main attributes are its very low cost per grid node and 
lack of dissipative error. Yee’s method is often applied using Cartesian grids, with a special 
tr 'itment of curved boundaries [5]. Extension to curvilinear grids was carried out by Fusco 
f Madsen and Ziolowski [7] put the method into a finite-volume framework applicable 
t‘ unstructured grids. Vinokur and Yarrow [8, 9] developed a related finite-surface method 
with advantages at boundaries and grid singularities. 
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Other methods which have been successfully applied to the time-domain Maxwell equa- 
tions include the upwind Lax-Wendroff approach used by Shankar [10], the characteristic- 
based fractional-step method of Shang [11], and the finite-element method of Cangellaris 
et al. [12] Although all of these second-order methods have been used with a great deal 
of success, they are efficient only for geometries of moderate electrical size, on the order 
of twenty wavelengths or less. For wave propagation over longer distances, the grid reso- 
lution requirements of second-order methods can become excessive, leading to impractical 
CPU and memory requirements. This has motivated the development of higher-order meth- 
ods which produce smaller errors for a given grid resolution, such as the extension of Yee’s 
method to fourth-order in space [2] and the methods of Liu [13], Shang [14], and Zingg et 
al. [15, 16, 17, 18] 

In seismology, the need for higher-order methods has been recognized for some time. 
Alford et al. [19], Marfurt [20], and Dablain [21] present higher-order algorithms for the 
elastic wave equation. Similarly, higher-order finite- difference methods have been developed 
for acoustic applications by Gottlieb and Turkel [22], Cohen and Joly [23], and Davis [24]. It 
can be advantageous to modify the coefficients of a potentially higher-order method, thereby 
lowering the order of accuracy, to produce reduced errors over a range of wavenumbers. This 
approach was first proposed by Vichnevetsky and De Schutter [25] and later studied in more 
detail by Holberg [26] and Lele [27], Haras and Ta’asan [28] and later Kim and Lee [29] 
further refined the optimization technique used by Lele. The papers by Holberg and Lele 

spawned a number of optimized schemes, including those presented in Refs [151 [161 and 
[30] to [34], 1 J 1 J 

Numerical errors arise from both the spatial and the temporal discretization. They 
include both phase and amplitude errors, which depend on the wavenumber, the grid spacing, 
the Courant number, and the direction of propagation relative to the grid. The dependence of 
the phase speed on the wavenumber results in numerical dispersion. This paper presents the 
grid resolution required to achieve a specific level of accuracy as a function of the propagation 
distance expressed in terms of the wavelength for several high-order and optimized methods. 
Emphasis is on methods requiring under 30 grid points per wavelength {PPW) for accurate 
simulations with propagation distances of 200 wavelengths. The errors are calculated using 
Fourier analysis [35, 36, 16]. The purpose is to aid in selecting and using a scheme for a 
given application, and to provide direction for future research. 

In addition to the accuracy of the interior differencing scheme, there are several other 
important issues in simulating wave phenomena. High-order and optimized methods often 
have a large spatial stencil which cannot be used at boundaries. This necessitates the use 
of numerical boundary schemes (NBS’s) which must be appropriately accurate and stable. 
These can be difficult to obtain, and thus this represents the most significant obstacle in the 
use of higher-order methods. Recent progress is reported in Refs. [37] to [40]. Another im- 
portant consideration is the boundary condition at the outer boundary of the domain, which 
inevitably causes spurious reflection. Since this does not decrease with the grid spacing, 
it can become the limiting factor in determining the accuracy of a simulation. Important 
early contributions to the development of non-reflecting boundary conditions are given by 
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Engquist and Majda [41] and Bayliss and Turkel [42]. Mur [43] implemented the Engquist- 
Majda approach for the time-domain Maxwell equations. Moore et al. [44] present a detailed 
comparison of non-reflecting boundary conditions available at that time, again in an elec- 
tromagnetic context. A review of developments up to 1991 is given by Givoli [45]. Progress 
since that time includes the approaches of Tam and Webb [46] and Berenger [47], both of 
which show considerable promise. A more recent comparison, in an aeroacoustic context, is 
presented by Hixon et al. [48] 

Another important consideration is the choice of a gridding strategy. Although un- 
structured grids have been used [7], Cartesian and curvilinear grids are generally preferred. 
Uniform Cartesian grids are used, for example, by Taflove [2] and Tam [49], while body-fitted 
curvilinear grids are favoured by Shankar et al. [10], Jurgens and Zingg [17, 18], and others. 
The Cartesian approach produces two significant advantages in the interior of the domain. 
Virtually all numerical methods are most accurate on a uniform grid, and furthermore, the 
grid size can be chosen to give the Courant number which produces the smallest errors for 
the scheme used. In addition, the need to deal with the metrics of a curvilinear coordinate 
transformation leads to a substantial penalty in terms of both speed and memory. On the 
other hand, it is extremely difficult to develop stable and accurate boundary treatments for 
higher-order methods on Cartesian grids. Refs. [5] and [50] present progress in this area. 
Furthermore, some geometries, such as a curved surface with a thin coating, are not well 
suited to Cartesian grids. Clearly, the choice of a gridding strategy is problem dependent. 


Finite-Difference Methods Compared 


In this section, we present the various finite-difference methods in the context of the linear 
advection equation given by 


du du 
dt dx 


= 0 


( 1 ) 


where u is a scalar quantity propagating with speed a, which is real and positive. The spatial 
difference operators are thus approximations to d/dx. Extension to systems of equations and 
multidimensions is not discussed. The time-marching methods are presented as applied to a 
scalar ordinary differential equation of the form 


du 

¥ = /K<) (2) 

The fully-discrete schemes are given as applied to eq. 1 itself. The coefficients of all of the 
schemes studied are given in the Appendix. 
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Spatial Difference Operators 

On a uniform grid with x 3 = jAx and uj = u(x.j), compact centered-difference schemes of 
up to tenth order can be represented by the following formula: 

0{bxu ) 3 - 2 + a(^ r u)j_j + (S x u)j + a(<^ x u)j +1 + 0(8 x u) 3+2 

1 c b 

= 2Az 3^ i+3 ~~ + 2^ Uj+2 ~ Uj ~ 2 ^ + a ^ j+1 ~ ( 3 ) 

where is an approximation to djdx at node j. These schemes produce no dissipative 

error. Noncompact schemes of up to sixth order are obtained with (3 = 8 = 0. 

Haras and Ta’asan [28] revisited the compact schemes of Lele [27] based on eq. 3 using 
a more sophisticated optimization procedure. They developed tridiagonal schemes (0 = 
0) and tridiagonal schemes with a five-point stencil (0 = c = 0) as well. In each case, 
Haras and Ta’asan present several methods which result from different parameters in the 
optimization procedure. In the comparisons below, we will include only those schemes which 
are best according to our criterion, i.e., those schemes which require the fewest grid points 
per wavelength for accurate wave propagation over a distance of two hundred wavelengths. 

The spatial operator of Zingg et al. [15, 16] is noncompact with a seven-point stencil. 
The operator is divided into an antisymmetric part given by 

(*»; = [ fl 3 ( u j+3 — Uj-3 ) + a 2 (u J+ 2 — Uj_ 2 ) +ai {Uj+i — Uj — i ) ] (4) 

and a symmetric part given by 

(fix u )j = [^3 (uj+3 + U J~3) + d 2 {u J+ 2 + nj- 2 ) +dj (u J+1 + Uj_x) + d 0 Uj\ (5) 

The symmetric part provides dissipation of spurious high- wavenumber components of the 
solution and is required for stability of the time-marching method used, which will be pre- 
sented in the next subsection. The magnitude of the dissipative component is chosen such 
that the dissipative error produced is comparable to the dispersive error. Refs. [15] and 
[16] include a maximum-order scheme and an optimized scheme. The optimized scheme was 
designed to minimize the maximum phase and amplitude errors for waves resolved with at 
least ten PPW. It is superior for distances of travel of up to 330 wavelengths. Hence we 
do not consider the maximum-order scheme here. Tam [31] has also developed an optimized 
scheme based on the antisymmetric operator given in eq. 4. 

Lockard et al. [30] present an optimized upwind-biased noncompact spatial difference 
operator based on an eight-point stencil. It can be written in the form 

1 3 

(8 x u)j = — — a m Uj+ m ( 6 ) 

m =— 4 

The operator is optimized for waves resolved with at least seven points per wavelength and 
provides excellent accuracy up to about 340 wavelengths of travel. 
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Differencing schemes on staggered grids are well suited to problems in which the time 
derivative of one variable depends on the spatial derivative of the other, and vice-versa. This 
is the case for the time-domain Maxwell equations or the Euler equations linearized about a 
reference state with zero velocity, for example. Centered staggered differencing schemes of 
up to sixth order can be written in the form 

( 8 x u)j = — [&l (uj + l/2 - T7-I/2) + b 2 (u j+ 3/2 - Uj-3/2) + b 3 (t/,+5/2 - Uj_ 5/2) ] ( 7 ) 

Time-Marching Methods 

There are many considerations involved in selecting a time-marching method, including 
efficiency, i.e., accuracy per unit computational effort, stability, and memory use. Since wave 
propagation problems are generally not stiff, explicit methods are appropriate. Thus the 
Adams-Bashforth and Runge-Kutta families are suitable candidates. Since a large portion 
of the computational effort is generally associated with evaluation of the derivative function, 
one can approximately assess efficiency by accounting for the number of derivative function 
evaluations per time step. Runge-Kutta methods require one derivative function evaluation 
per stage. Zingg and Chisholm [51] have shown that for linear ordinary differential equations 
(ODE’s) with constant coefficients, Runge-Kutta methods of up to sixth-order (and probably 
arbitrary order) can be derived with a number of stages equal to the order. However, the 
memory requirements of the methods of order five and six are high. An alternative approach 
is to use low-storage multi-stage methods which are high-order for homogeneous linear ODE’s 
but second-order otherwise [52]. Haras and Ta’asan [28] and Zingg et al. [15, 16] present 
five- and six-stage methods of this type, respectively. For example, when applied to eq. 2, 
the six-stage method of Refs. [15] and [16] is given by 

= u n + ha\f n 
= u„ + ha 2 fn+ ai 

= u n + ha 3 f$ Q2 (8) 

= u n + ha 4 fn + a3 

= u n T /io;5/)1+q 4 
= u n + hf$ Qs 

where h = At is the time step, t n = nh, u n = u(t n ), and 

fn+ a = f{Unl a ,tn + ah) 

The five-stage method of Ref. [28] is analogous. In their maximum-order form (for homo- 
geneous ODE’s), these methods are not stable for pure central differencing in space. Con- 
sequently, Haras and Ta’asan modified the coefficients of the five-stage scheme to produce 
stability, while Zingg et al. added some dissipation to their spatial operator. 


u (1) 

U 7l + C*l 
( 2 ) 

u n+a 2 

( 3 ) 

n+a 3 

( 4 ) 

n+or 4 

(5) 

n+as 


U 


U 


U 


^n + 1 
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Since Adams-Bashforth methods require only one derivative function evaluation per time 
step, they are more efficient than Runge-Kutta methods of equal order. Unfortunately, 
Adams-Bashforth methods of order higher than four have extremely restrictive stability 
bounds which outweigh this advantage. Thus the time-marching methods we consider here 
are the fourth-order Adams-Bashforth method, the fourth-order Runge-Kutta method for lin- 
ear ODE’s given by Zingg and Chisholm [51], and the low-storage five- and six-stage methods 
described above. In terms of stability and Fourier error analysis, the fourth-order Runge- 
Kutta method of Ref. [51] is identical to the classical fourth-order Runge-Kutta method. 
With respect to memory requirements, the fourth-order Adams-Bashforth method requires 
five memory locations per dependent variable. The other methods considered require only 
two. 

A natural choice of time-marching method for use with staggered spatial differencing is 
the second-order staggered leapfrog method, given by 

u n+l = U n + /j/n+i/2 (9) 

When used with a nondissipative spatial scheme, as in the FDTD method, the resulting fully- 
discrete operator produces no dissipative error. Furthermore, the staggered leapfrog method 
generally produces a leading phase error, which can offset the phase lag usually produced 
by centered spatial differences. At specific Courant numbers and angles of propagation, 
the perfect-shift property can be obtained, leading to exact propagation for all wavenum- 
bers. Although this has little practical significance, fully-discrete methods which possess this 
property are generally most accurate when used near the perfect-shift conditions. 


Simultaneous Space-Time Discretizations 


The schemes presented in this subsection approximate the spatial and temporal derivatives 
simultaneously. No intermediate semi-discrete form is produced. All of the schemes are 
fourth-order in time and space. 

The scheme of Davis [24] is a nondissipative implicit scheme with a three-point spatial 
stencil. When applied to the linear advection equation, it is given by 


a 0 w" +1 + ain"^ 1 + a 2 u]+l = b 0 u n 3 + + b 2 u] +1 (10) 

with ao — bo , aj = b 2 , and a 2 = b\, where u” = t n ). Being implicit, this scheme requires 
somewhat more computational expense than explicit noncompact schemes. However, it 
achieves fourth-order accuracy with a three-point spatial stencil. 

The scheme of Gottlieb and Turkel [22] is an extension of the Lax-Wendroff approach 
to fourth-order. It is an explicit one-step scheme but requires a five-point stencil. When 
applied to the linear advection equation, it is given by 


u 


(i) _ 


ah 


u 


3 

n+ 1 


- + 6Ax^ +2 ^ u 


= iW + vW + ^-W + Suyi-ru?’) 


ah 


(i) 


, 0 ) 




12Ax 


( 11 ) 
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This method is dissipative and has received considerable use in nonlinear problems. 

Finally, we consider the following two-step explicit scheme [53]: 

' = u?- 1 - 2o*(4‘>)J - ^p(«£L)J (12) 

where 6^ is a fourth-order centered difference approximation to a first derivative and b^xx 
is a second-order centered approximation to a third derivative. Both of the operators on the 
right-hand side require a five-point stencil, as follows: 

( 4 4) )t = J2K1 ( Uj ~ 2 ~ + 8uj+1 ~ Uj+ ^ 

(^rrx)j = 2 Ax U j~ 2 ^ U j ~ 1 — ^ u j + 1 T 2 ) 

This scheme is derived using simple Taylor series expansions and is nondissipative. Related 
methods for the second-order wave equation are presented by Cohen and Joly [231 and Shubin 
and Bell [54], 


Fourier Error Analysis 

Fourier analysis provides a simple means of determining the errors produced by finite- 
difference methods in the absence of boundary conditions. In one dimension, these are 
functions of the wavenumber, k, the grid spacing, and the Courant number, C = ah/ Ax. 
In multidimensions, the errors further depend on the direction of propagation relative to 
the grid. The error analysis is based on the principal root, a(z, C), where z = kAx, i.e., 
z = 2i t/PPW. The spurious roots are required for stability analysis. The local amplitude 
and phase errors are, respectively, 


er a = |cr| — 1 (13) 

= (14) 

where <f> = arctan(<7{/<r r ) and a r and cr, are the real and imaginary parts of a, respectively. 
Our criterion for comparing schemes is based on the global amplitude and phase errors, 
which are 


Er 


a 


i _ PPWn„/C 


(15) 


Er p — ^hu 


^ + 2 , 


( 16 ) 
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where n w is the number of wavelengths travelled. Using the above formulas with a very 
small Courant number gives the errors for the spatial operator alone. In the figures below, 
the various methods are compared in terms of the PPW required to keep both global phase 
and amplitude errors below 0.1 as a function of the number of wavelengths travelled. 

The error from a spatial discretization is often plotted in terms of the modified wavenum- 
ber. Plots of the local phase and amplitude errors can be much more revealing and provide a 
stronger physical connection. The global errors are even easier to interpret. In studying the 
dependence of the PPW requirements on the number of wavelengths travelled, the emphasis 
is on the wavenumber present in the simulation which is most poorly resolved among those 
wavenumbers which are deemed to be significant. This is a reasonable measure for selecting 
a grid density. Furthermore, this method of presentation clearly reveals the implications of 
optimization. 

Time-marching methods can be analyzed in isolation by applying them to the following 
linear homogeneous ODE: 

du 

li = Xu < 17 > 

In order to assess their accuracy for wave propagation, we consider only pure imaginary 
values of A, i.e., A = ito with 10 real. The local amplitude and phase errors are determined 
from the principal root, <t(A/j), as follows: 


er a = \c\ — 1 


(18) 


er„ 


arctan(<T,-/<7 r ) 

ujh 


+ 1 


(19) 


Results 


Spatial Difference Operators 

In this section, we consider the errors produced by the spatial operators alone. Figure 1 shows 
the grid resolution requirements for second-, fourth-, and sixth-order noncompact centered 
difference schemes. Since these spatial operators are nondissipative, the grid resolution 
requirements are determined from the global phase error. The requirements of the second- 
order scheme are clearly excessive. While more accurate second-order discretizations are 
available, such as the FDTD scheme or the upwind leapfrog scheme [55], none of these meets 
our criterion of 30 PPW for 200 wavelengths of travel except under specific conditions [56]. 

The PPW requirements for compact centered difference schemes of up to tenth order are 
shown in Figure 2. The tenth-order scheme requires the solution of a pentadiagonal system 
of equations. The remaining schemes lead to tridiagonal systems. The compact schemes are 
considerably more accurate than the noncompact schemes of equivalent order. 
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Figure 1: Grid resolution requirements for second-order ( ), fourth-order ( ), and sixth- 

order (• ■ •) noncompact centered spatial differences. 



Figure 2: Grid resolution requirements for fourth-order ( ), sixth-order ( ), eighth-order 

(- • •), and tenth-order (— • — ) compact centered spatial differences. 
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Figure 3: Grid resolution requirements for the tridiagonal spatial operators with 5-point 

right-hand stencil of Haras and Ta’asan: scheme A ( ), scheme B ( ), scheme C (• • •), 

sixth-order scheme (— • — ). 


We next demonstrate the tradeoffs associated with optimization, using the tridiagonal 
five-point operators (/? = c = 0) of Haras and Ta’asan as an example. Figure 3 shows the 
behavior of three different optimized schemes given in Ref. [28] as well as the maximum- 
order scheme which can be obtained using this operator, which is sixth-order. The scheme 
denoted “A” is superior up to a distance of travel of about 60 wavelengths. Scheme B is 
superior for distances up to 250 wavelengths while scheme C is preferable for even longer 
distances. Such behavior is typical of optimized schemes. Aggressive optimization (as in 
the case of scheme A) leads to excellent performance for small distances of travel but poor 
performance for longer distances. In this paper we concentrate on distances of travel up to 
200 wavelengths. Hence we consider only scheme B further. 

Figure 4 shows the PPW required by the optimized schemes of Haras and Ta’asan for 
distances of travel up to two hundred wavelengths. In each case, the best scheme presented 
by Haras and Ta’asan for this distance of travel is shown, as discussed above. The penta- 
diagonal scheme requires about 3.7 PPW for two hundred wavelengths of travel while the 
tridiagonal seven-point scheme requires 4.5 and the tridiagonal five-point scheme requires 
5.7. The computational effort is roughly proportional to PPW d+i , where d, is the number 
of dimensions. Therefore the increased cost per grid node of the pentadiagonal scheme is 
probably not justified. However, the tridiagonal seven-point scheme is likely to be more 
efficient than the tridiagonal five-point scheme. 

The noncompact optimized schemes of Lockard et ah, Zingg et ah, and Tam are compared 
with sixth-order centered differences in Figure 5. For the scheme of Lockard et ah, the PPW 
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Figure 4: Grid resolution requirements for the spatial operators of Haras and Ta’asan: pen- 

tadiagonal operator with 7-point right-hand stencil ( ), tridiagonal operator with 7-point 

right-hand stencil ( ), tridiagonal operator with 5-point right-hand stencil (• • •). 


requirements are determined by the amplitude error. For the scheme of Zingg et al., the 
phase and amplitude errors produce roughly equivalent grid resolution requirements over 
this distance. Since Tam’s scheme is nondissipative, its PPW requirements are determined 
by the phase error. The figure shows that Tam’s scheme has been optimized for fairly small 
distances of travel, leading to excellent performance for less than roughly 14 wavelengths of 
travel. The scheme of Lockard et al. requires roughly 8 PPW for two hundred wavelengths 
of travel while that of Zingg et al. requires over 9 PPW . This again is probably sufficient 
to compensate for the increased expense of the Lockard scheme, which uses an eight-point 
stencil versus seven for the Zingg scheme. However, it can be difficult to find stable and 
accurate numerical boundary schemes for use near boundaries when the spatial stencil is 
large. For the scheme of Zingg et al., such numerical boundary schemes are given in Ref. [17]. 
Comparison of these two schemes with those of Haras and Ta’asan is difficult. Although the 
compact schemes have significantly reduced grid resolution requirements, they also require 
considerably greater computational effort per grid node. Numerical experiments are required 
to determine which approach is more efficient. 

Figure 6 shows the PPW requirements of the staggered spatial operators. In each case, 
the staggered schemes are much more accurate than their nonstaggered counterparts of 
equivalent order. The grid requirements of the second-order scheme are again excessive. 
However, when used with staggered leapfrog time marching (the FDTD scheme), better 
results can be obtained. Also shown in Figure 6 is an optimized scheme with i 3 = 103/19200, 
f >2 = —1315/19200, and b\ = 22630/19200, which produces excellent performance for up to 
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number of wavelengths travelled 

Figure 5: Grid resolution requirements for the spatial operators of Lockard et al. ( ), 

Zingg et al. ( ), Tam (■ • •), and sixth-order centered differences (—•—). 

200 wavelengths of travel. 


Time-Marching Methods 

Figures 7 and 8 show the amplitude and phase errors produced by the four time-marching 
methods under consideration. The five- and six-stage methods shown are the maximum-order 
versions rather than the optimized methods, which are discussed in the next subsection. In 
order to account for the number of stages in the Runge-Kutta methods, the errors are plotted 
versus uh/p, where p is the number of stages. Hence the errors shown are for approximately 
equal computational effort. Since the time step of a p-stage scheme is thus p times larger 
than that of a single-stage scheme, the amplitude error shown is |<r| 1 / p — 1. 

The figures show that the phase errors produced by these methods are larger than the 
amplitude errors. Each increase in the order of the Runge-Kutta method produces an increase 
in accuracy even though the extra work has been accounted for. The fourth-order Adams- 
Bashforth method is much more accurate than the fourth-order Runge-Kutta method. It 
produces the lowest amplitude error of the methods considered and phase error comparable 
to the five-stage Runge-Kutta method. 

In order to compare time-marching methods properly, one must consider the spatial 
discretization to be used. For each combination of a spatial discretization and a time- 
marching method, there is a Courant number which minimizes the computational work to 
achieve a given standard of accuracy for a specified distance of travel. The computational 
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number of wavelengths of travel 


Figure 6: Grid resolution requirements for second-order ( ), fourth-order ( ), sixth-order 

(• • •), and optimized (— • — ) centered spatial differences on a staggered grid. 



(o h/p 

Figure 7: Amplitude error produced by the fourth-order Adams-Bashforth method ( ), 

the fourth-order Runge-Kutta method ( ), the five-stage Runge-Kutta method (• • •), and 

the six-stage Runge-Kutta method (— • — ). 
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Figure 8: Phase error produced by the fourth-order Adams-Bashforth method ( ), the 

fourth-order Runge-Kutta method ( ), the five-stage Runge-Kutta method (• • ■), and the 

six-stage Runge-Kutta method (— • — ). 


work is proportional to 

F = £ PPW d+1 (20) 

where p is again the number of stages and d the number of dimensions. For example, 
with fourth-order centered differences in space the optimum Courant number for the fourth- 
order Runge Kutta method in two dimensions is about 1.25 for two hundred wavelengths 
of travel and global errors less than 0.1. With this spatial operator, the optimum Courant 
number for the Adams-Bashforth method is determined by its stability bound, which for 
pure imaginary A is roughly 0.43. With fourth-order centered differences in two dimensions, 
the resulting Courant number is about 0.21. The value of F produced at this Courant 
number is slightly higher than that produced using the fourth-order Runge-Kutta method at 
its optimum Courant number. Since the fourth-order Runge-Kutta method also requires less 
memory than the fourth-order Adams-Bashforth method, it is clearly the preferred method 
for use with fourth-order centered differences in space. 

With sixth-order centered differences in space, the comparison changes. The optimum 
Courant number for the fourth-order Runge-Kutta method is reduced to roughly 2/3 while 
that for the fourth-order Adams-Bashforth method remains stability limited, leading to a 
Courant number of roughly 0.19. In this case the value of F produced by the Adams- 
Bashforth method is less than 80% of that for the Runge-Kutta method. Thus the compari- 
son is more complicated, as one must weigh the increased efficiency of the Adams-Bashforth 
method against the reduced memory requirements of the Runge-Kutta method. 

The difficulty with the maximum-order five- and six-stage Runge-Kutta methods is that 
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they are unstable for pure imaginary A, as shown in Figure 7. We consider these schemes 
further below. 

Combined Space-Time Discretizations 

Haras and Ta’asan modified the coefficients of the five-stage Runge-Kutta method to produce 
stability for pure imaginary A while maintaining second-order accuracy and optimized error 
behavior. The methods are designed for C = 0.9. Figure 9 shows the grid resolution 
requirements of the three spatial operators of Haras and Ta’asan compared in Figure 4 in 
combination with their five-stage time-marching method at a Courant number of 0.9. All 
three schemes require between 7 and 8 PPW for two hundred wavelengths of travel. The 
advantage of the more accurate spatial operators has been lost. Either a lower Courant 
number or a more accurate time-marching method should be used. 

As a result of the dissipation in the spatial operator, the six-stage time-marching method 
of Zingg et al. is stable up to a Courant number a little greater than unity in two dimensions. 
The grid requirements for the combined space-time discretization at a Courant number of 
unity are shown in Figure 9. Comparison with Figure 5 shows that the time-marching 
method introduces very little error compared to the spatial differencing. Tam [31] uses an 
optimized four-step Adams-Bashforth method in conjunction with his spatial operator, ft 
produces little error for Courant numbers less than 0.3. In both cases, optimization of the 
time-marching method has a much smaller impact than optimization of the spatial operator. 

Figure 10 shows the PPW requirements of the fourth-order staggered spatial difference 
operator combined with staggered leapfrog time marching. As the Courant number is in- 
creased from 0.001 to 0.1, the PPW requirements decrease, since the error from the time 
marching method is of opposite sign to that of the spatial operator. However, as the Courant 
number is further increased, the error from the second-order time-marching method begins 
to dominate. Excellent performance for 200 wavelengths of travel is obtained for Courant 
numbers up to 0.1. 

Results for the method of Davis are shown in Figure 11. Since the method is nondissi- 
pative, the PPW requirements are determined from the phase error. This scheme is uncon- 
ditionally stable. For Courant numbers under 1.5, less than 19 PPW are required for two 
hundred wavelengths of travel. While this is quite good, much better than many schemes, 
it is probably not sufficient to justify the additional computational effort associated with an 
implicit scheme. 

Figure 12 shows the grid resolution requirements for the Gottlieb-Turkel scheme, which is 
stable up to a Courant number of 2/3 in one dimension. This scheme is very robust and has 
received extensive use in nonlinear applications. However, high accuracy is obtained only for 
Courant numbers less than about 0.13. Less than 29 PPW are required to propagate a wave 
two hundred wavelengths with Courant numbers below this value. Although the cost per 
grid node is reasonably low, these PPW requirements are too high for efficient simulation 
over such distances. 

Finally, the grid requirements of the nondissipative two-step explicit scheme, eq. 12, 
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Figure 9: Grid resolution requirements for the spatial and temporal operators of Haras and 
Ta’asan at a Courant number of 0.9: pentadiagonal operator with 7-point right-hand stencil 

( )> tridiagonal operator with 7-point right-hand stencil ( ), tridiagonal operator with 

5-point right-hand stencil (• • •); spatial and temporal operators of Zingg et al. at a Courant 
number of 1 (— • — ). 



Figure 10: Grid resolution requirements for fourth-order centered differences on a staggered 

grid coupled with staggered leapfrog time marching at Courant numbers of 0.2 ( -), 0.1 (- 

- -), and 0.001 (• • •). 
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number of wavelengths of travel 


Figure 11: Grid resolution requirements for the scheme of Davis at Courant numbers of 3 
( ), 1.5 ( ), and 0.01 (• ■ •). 



number of wavelengths of travel 


Figure 12: Grid resolution requirements for the scheme of Gottlieb and Turkel at Courant 
numbers of 0.2 ( ), 0.13 ( ), and 0.01 (■ • •). 
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Figure 13: Grid resolution requirements for the two-step explicit scheme (eq. 10) at Courant 
numbers of 0.9 ( ), 0.5 ( ), and 0.01 (•• •). 

are shown in Figure 13. This scheme is stable up to a Courant number of unity in one 
dimension. It has the perfect shift property at a Courant number of unity. The phase error 
increases as the Courant number is reduced. Less than 29 PPW are required for all stable 
Courant numbers. Since it provides reasonably low error at Courant numbers up to unity, 
this scheme has a very low cost per grid node, substantially lower than the Gottlieb-Turkel 
scheme. However, the PPW requirements are again much higher than some of the other 
schemes considered. 


Discussion 

The grid requirements presented can be used to determine a suitable grid density and Courant 
number when applying a method to a specific problem. Estimates of the shortest wavelength 
of interest and the distance of travel are required. However, long distances of travel are 
generally associated with multiple reflections. Therefore, the accuracy of a simulation is 
strongly dependent on the accuracy of the boundary condition formulation as well as the 
interior finite-difference scheme. This represents a major obstacle in the application of high- 
accuracy schemes such as those of Haras and Ta’asan. 

The results show that optimized schemes provide a significant advantage over their 
maximum-order counterparts. However, if optimized too aggressively, they perform poorly 
for longer distances of travel. If a waveform has significant low wavenumber content, as 
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in the case of a Gaussian pulse, maximum-order schemes can produce smaller errors than 
optimized schemes. Furthermore, numerical boundary schemes affect the relative accuracy 
of optimized and maximum-order schemes. 

The spatial operators of Haras and Ta’asan, Lockard et ah, and Zingg et al. all provide 
adequate accuracy for two hundred wavelengths of travel. The selection of a time-marching 
method is not quite as critical, since the computational work varies linearly with the time 
step size. Consequently, the cost of reducing the time step is much less than the cost of 
reducing the PPW , and, furthermore, has no memory implications. Adams-Bashforth and 
low-storage Runge-Kutta methods can be used, with the latter often preferable because of 
their reduced memory requirements. Optimization of the time-marching method is generally 
not particularly beneficial. 

For appropriate problems, such as those involving electromagnetic waves or acoustic 
waves in a quiescent medium, staggered spatial schemes perform very well. These can be 
combined with either a high-order time-marching method or the staggered leapfrog method. 
In the latter case, a low Courant number should be used. For large propagation distances, 
higher-order time-marching methods are more efficient. 

The grid requirements of the explicit fourth-order methods involving simultaneous space- 
time discretization are reasonably low considering the low cost of these schemes, especially 
the two-step explicit scheme. This suggests that sixth-order extensions of these schemes are 
worthy of investigation. Optimized forms of these schemes can be considered as well. 


A Numerical Example 

In this section, we present a numerical simulation of the propagation and reflection of an 
electromagnetic wave which demonstrates some of the issues discussed. Further details are 
available in Ref. [18]. The governing equations are the transverse magnetic set of the two- 
dimensional time-domain Maxwell equations. The simulation consists of a pulsed plane 
wave incident on a perfectly-conducting cylinder. A grid containing 5,400 nodes is shown 
in Figure 14. Figure 15 shows a snapshot of the electric field intensity computed on a grid 
with 21,600 nodes using the maximum-order version of the method of Zingg et al. [15, 16] 
The dashed contours indicate negative values of the electric field intensity. This solution 
is visually indistinguishable from that computed using the same method on a grid with 16 
times as many nodes, which is used as a baseline to estimate numerical errors. 

Figures. 16 and 17 show the Li norm of the numerical error versus the number of 
grid nodes and the CPU effort, respectively. Four methods are included, the maximum- 
order (MO) and optimized (OlO) schemes of Zingg et al. [15, 16], as well as second- (C2) 
and fourth-order (C4) centered differences. The second- and fourth-difference schemes are 
combined with fourth-order Runge-Kutta time marching and include a very small amount 
of artificial dissipation for stability. These figures show that the higher-order and optimized 
methods lead to substantial reductions in both memory and CPU time. The benefits of the 
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Figure 15: Computed contours of electric field intensity. 
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optimized scheme over the maximum-order scheme are fairly modest due to the nature of 
the waveform, which is a Gaussian and thus has considerable low-wavenumber content. 


Conclusions 

Based on the results presented, it is clear that high-order and optimized finite-difference 
methods will play an important role in the simulation of high frequency linear wave phe- 
nomena. Several of the methods studied have the potential to reduce substantially the 
computational requirements, including both CPU time and memory, for accurate simula- 
tions. Further study is required to determine the relative advantages of these methods in a 
practical context. Development of stable and accurate numerical boundary schemes remains 
a pacing item in the application of such methods. Further effort is also required in the devel- 
opment of gridding strategies permitting high-accuracy methods to be applied to complex 
geometries. 
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Appendix 

The following are the coefficients of the finite-difference schemes studied. 

The spatial operator of Haras and Ta’asan is given in eq. 3. The schemes shown in Figure 
3 have ft = c = 0. The remaining coefficients are: 

Scheme A: 

a = 0.3534620453 
a = 1.566965775 
b = 0.1399583152 

Scheme B: 

a = 0.3461890571 
a = 1.5633098070 
b = 0.1290683071 
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Scheme C: 


a = 0.3427812069 
a = 1.5614141543 
b = 0.124148259 

In Figure 4, the pentadiagonal 7-point scheme has 

a = 0.5801818925 
P = 0.0877284887 
a = 1.3058941939 
b = 0.9975884963 
c = 0.0323380724 

The tridiagonal 7-point scheme has 

a = 0.3904091387 

P = 0 

a = 1.5638887738 
b = 0.2348222711 
c = -0.0178927675 

The tridiagonal 5-point operator is scheme B above. 

In Figure 5, the scheme of Lockard et al. is the average of two schemes with the following 
coefficients, as defined in eq. 6: 


and 


a_ 4 = 0.0207860419 
G_ 3 = -0.1500704734 
0—2 = 0.5234309723 
0—i = -1.34207332539 
G 0 = 0.574548248808 
ai = 0.4090357053658 
a 2 = -0.035657169508 


a Q = 

0 


a\ — 

-o-i = 

0.763289242273 

a 2 = 

—0-2 = 

-0.160631393818 

«3 = 

-0-3 = 

0.019324515121 



The optimized scheme of Zingg et al. is given by eqs. 4 and 5 with the following coefficients: 

oi = 0.75996126 
a 2 = -0.15812197 
a 3 = 0.018760895 
d 0 = 0.1 

di = -0.076384622 
d 2 = 0.032289620 
d 3 = -0.0059049989 

Tam’s scheme is obtained from eq. 4 with 


a t = 0.770882380518 
a 2 = -0.166705904415 
a 3 = 0.0208431427703 

The five-stage Runge-Kutta method shown in Figures. 7 and 8 has the following charac- 
teristic polynomial: 
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The six-stage method is obtained from eq. 8 with a 5 = 1/2, a 4 = 1/3, a 3 = 1/4, a 2 = 1/5, 
aj = 1/6, leading to the following characteristic polynomial: 
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The characteristic polynomial of the optimized five-stage Runge-Kutta method of Haras 
and Ta’asan used in Figure 9 is: 


a{Xh) = 1 4- Xh + + 0.166407(A/>) 3 

+0.0409525(A/i) 4 + 0.0074510(A/t) 5 

The optimized temporal operator of Zingg et al. is obtained from eq. 8 with the following 
coefficients: 


a] = 0.168850 
q 2 = 0.197348 
q 3 = 0.250038 
a 4 = 0.333306 
a 5 = 0.5 
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The resulting characteristic polynomial is 

a(Xh) = 1 + A h + + 0.16665295(A/i) 3 

+0.041669557(A/i) 4 + 0.0082233848( A*) 5 
+0.0013885169 (A/j) 6 

The scheme of Davis is obtained from eq. 10 with 

Go = b 0 = — 2(C — 2)(C + 2) 

fli = b 2 = (C — 1)(C — 2) 

a 2 — bi = ( C + 1)(C -f 2) 

where C = ah/ Ax is the Courant number. 


29 




